#version 420

uniform sampler3D TexDiagTensor;
uniform sampler3D TexOffDiagTensor;

in vec3 pos;
in vec3 norm;
in vec3 texCoord;

out vec4 frontColor;

struct Light {
    vec4 pos;        // Light position in eye coords.
    vec3 ambient;    // Ambient light intensity
    vec3 diffuse;    // Diffuse light intensity
    vec3 specular;   // Specular light intensity
};

struct Material {
    vec3 ka;         // Ambient reflectivity
    vec3 kd;         // Diffuse reflectivity
    vec3 ks;         // Specular reflectivity
    float shininess; // Specular shininess factor
};

struct Vertex {
    vec3 pos;
    vec3 normal;
    Material material;
};

Light Glight = Light(
    vec4(-1.0, -1.0, -1.0, 0.0),
    vec3(1.0, 1.0, 1.0),
    vec3(1.0, 1.0, 1.0),
    vec3(1.0, 1.0, 1.0)
    );

vec3 phong(in vec3 antiViewDirection/*opposite view direction*/, in Light light, in Vertex vertex)
{
    vec3 antiLightDirection;    // opposite incident light direction


    if(0.0 == light.pos.w)  //  directional light
    {
        antiLightDirection = normalize(-vec3(light.pos));
    }
    else
    {
        antiLightDirection = normalize(light.pos.xyz - vertex.pos);
    }

    vec3 ambient = light.ambient * vertex.material.ka;
    vec3 diffuse = light.diffuse * vertex.material.kd * dot(antiLightDirection, vertex.normal);

    vec3 reflectLight = reflect(antiLightDirection, vertex.normal);
    float kReflect = pow(dot(reflectLight, antiViewDirection), vertex.material.shininess);
    
    //vec3 hv = antiLightDirection + antiViewDirection;
    //float kReflect = pow(dot(hv, vertex.normal), vertex.material.shininess);
    
    vec3 specular = light.specular * vertex.material.ks * kReflect;

    return clamp(ambient + diffuse + specular, vec3(0.0), vec3(1.0));
}

/*
 * \note:
 * tensors have to be stored in the following way:
 * diag.x = Dxx     Diagonal elements are stored in one vector
 * diag.y = Dyy
 * diag.z = Dzz
 * offdiag.x = Dyz  Off diagonal elements are stored in another vector
 * offdiag.y = Dxz      where all components are stored in that location
 * offdiag.z = Dxy      that is indexed by the coordinate that is not in the tensor index
*/

float sqr( float a )
{
    return a * a;
}

vec3 getEigenvaluesCardano( in vec3 diag, in vec3 offdiag )
{
  const float M_SQRT3 = 1.73205080756887729352744634151;
  float de = offdiag.z * offdiag.x;
  float dd = sqr( offdiag.z );
  float ee = sqr( offdiag.x );
  float ff = sqr( offdiag.y );
  float m  = diag.x + diag.y + diag.z;
  float c1 = diag.x * diag.y + diag.x * diag.z + diag.y * diag.z
             - ( dd + ee + ff );
  float c0 = diag.z * dd + diag.x * ee + diag.y * ff - diag.x * diag.y * diag.z - 2. * offdiag.y * de;
  float p, sqrt_p, q, c, s, phi;
  p = sqr( m ) - 3. * c1;
  q = m * ( p - ( 3. / 2. ) * c1 ) - ( 27. / 2. ) * c0;
  sqrt_p = sqrt( abs( p ) );

  phi = 27. * ( 0.25 * sqr( c1 ) * ( p - c1 ) + c0 * ( q + 27. / 4. * c0 ) );
  phi = ( 1. / 3. ) * atan( sqrt( abs( phi ) ), q );

  c = sqrt_p * cos( phi );
  s = ( 1. / M_SQRT3 ) * sqrt_p * sin( phi );

  vec3 w;
  w[2] = ( 1. / 3. ) * ( m - c );
  w[1] = w[2] + s;
  w[0] = w[2] + c;
  w[2] -= s;
  return w;
}

vec3 getEigenvector( vec3 ABC /*diag without eigenalue i*/, vec3 offdiag )
{
  vec3 vec;
  vec.x = ( offdiag.z * offdiag.x - ABC.y * offdiag.y ) * ( offdiag.y * offdiag.x - ABC.z * offdiag.z );
  vec.y = ( offdiag.y * offdiag.x - ABC.z * offdiag.z ) * ( offdiag.y * offdiag.z - ABC.x * offdiag.x );
  vec.z = ( offdiag.z * offdiag.x - ABC.y * offdiag.y ) * ( offdiag.y * offdiag.z - ABC.x * offdiag.x );

  return normalize( vec );
}

float getFA( vec3 evalues )
{
  float sqrt32 = sqrt( 3. / 2. );
  float I1 = evalues.x + evalues.y + evalues.z;
  float mu = I1 / 3.;
  vec3 deriv;
  deriv.x = ( evalues.x - mu );
  deriv.y = ( evalues.y - mu );
  deriv.z = ( evalues.z - mu );
  float FA = sqrt32 * length( deriv ) / length( evalues );
  return FA;
}

float getFA( vec3 diag, vec3 offdiag )
{
  float cross = dot( offdiag, offdiag );
  float j2 = diag.x * diag.y + diag.y * diag.z + diag.z * diag.x - cross;
  float j4 = dot( diag, diag ) + 2. * cross;
  return sqrt( ( j4 - j2 ) / j4 );
}


float calculateFA(in vec3 eigenValue)
{
    float RA = (eigenValue.x + eigenValue.y + eigenValue.z) / 3.0;

    vec3 d = eigenValue - vec3(RA);
    float u = d.x*d.x + d.y*d.y + d.z*d.z;

    return sqrt(1.5*u) / length(eigenValue);
}

// eigenVector should be normalized
vec3 colorEncoded(in vec3 value, in vec3 vector)
{
    float sum = value.x + value.y + value.z;

    // linear
    float cl = (value.x - value.y) / sum;
    // planar
    float cp = 2.0*(value.y - value.z) / sum;
    // spherical
    float cs = 3.0*value.z / sum;

    float u1 = sum / 3.0;
    vec3 dif = value - vec3(u1);

    vec3 tmp = pow(dif, vec3(2.0));
    float u2 = (tmp.x + tmp.y + tmp.z) / 3.0;

    tmp = pow(dif, vec3(3.0));
    float A3 = (tmp.x + tmp.y + tmp.z) / (3.0 * u2 * sqrt(u2));

    tmp = pow(value, vec3(2.0));
    float J4 = tmp.x + tmp.y + tmp.z;

    float FA = 3.0 * sqrt(u2 / J4 / 2.0);

    float RA = sqrt(u2 / 2.0) / u1;

    vec3 normVec = abs(vector);

    //vec3 RGB = vec3(u1, u2, A3);
    //vec3 RGB = vec3(J4, FA, A3);
    vec3 RGB = FA*normVec;

    return RGB;
}

void main()
{
    vec3 diagnal = texture(TexDiagTensor, texCoord).xyz;
    vec3 offDiag = texture(TexOffDiagTensor, texCoord).xyz;

    vec3 eigenValue = getEigenvaluesCardano(diagnal, offDiag);
    vec3 eigenVector = getEigenvector(diagnal-eigenValue.x, offDiag);

    vec3 color = colorEncoded(eigenValue, eigenVector);

    Vertex vertex = Vertex(
        pos,
        norm,
        Material(0.7*color, color, color, 20));

    color = phong(vec3(0.0, 0.0, 1.0), Glight, vertex);

    frontColor = vec4(color, 1.0);
}
